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ABSTRACT 

There is apparently a widespread belief that the gravitational field (and sub- 
sequently the rotation curve) "inside" razor-thin, axially symmetric disks can not 
be determined accurately from elliptic integrals because of the singular kernel in 
the Poisson integral. Here, we report a simple and powerful method to achieve 
this task numerically using the technique of "density splitting". 

Subject headings: gravitation — methods: numerical — Galaxy: disk — Galaxy: 
kinematics and dynamics 



1. Introduction. Statement of the problem 

One of the oldest problem in galactic dynamics is the determination of the mass dis- 
tribution in the stellar and gaseous disk from observed velocities. Observationally, noisy 
measurements along with a certain velocity dispersion makes this inverse problem difficult 
to solve (Beauvais & Bothum, 2001). From a theoretical point of view, various levels of 
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approximation are possible. For instance, the inversion is straightforward if the disk is con- 
sidered as a flattened spheroid, since the surface density is simply obtained by an inverse Abel 
transform over the velocity field (Brandt, 1960). Flattened spheroids are however a poor 
representation of real disks (e.g. Kochanek 2001). Another, more intuitive way to proceed 
is to generate a collection of surface density profiles, then to determine the corresponding 
potential or gravitational field, and find only those which model the velocity data at best. 
For a razor-thin axi-symmetrical disk, only one numerical quadrature over its radial extent 
must be performed to obtain the potential or the field, which is a priori very convenient be- 
cause time in-expansive. Meanwhile, this operation is commonly believed to be numerically 
uncomfortable when using elliptic integrals (Binney & Tremaine, 1987, Cuddeford, 1993, 
Kochanek, 2001). The reason is that these functions contain a singularity everywhere in the 
source. As suggested by Binney & Tremaine (1987), this drawback can be simply avoided by 
considering field points located just above/below the equatorial plane. Even if the potential 
and its radial derivative do generally not exhibit strong variations slightly off the mid-plane, 
such a vertical shift introduces an inconsistency or a bias. As a matter of fact, instead 
of elliptic integrals, most authors work with Bessel functions which have finite amplitudes 
(Toomre, 1963). But their oscillatory behavior and specially their spatial extension (the 
infinite range) are a severe disadvantage from a numerical point of view (Cuddeford, 1993; 
Cohl & Tohline 1999). 

In this paper, we show that point mass singularities can be properly and easily handled 
numerically from elliptic integrals by density splitting. For the case of a disk with zero 
thickness as considered here, the gravitational potential and field can be determined exactly 
in the equatorial plane without any shift, and whatever the surface density profile (provided 
it vanishes at the inner and outer edges of the disk). The extension to tri-dimensional disks 
is possible. 



For a disk with zero thickness, the expression for the radial gravitational field in the 
equatorial plane is known in a closed form. Using cylindrical coordinates, this reads (e.g 
Durand, 1964): 



where £ is the total surface density, r is the field point, a refers to the material distri- 
bution, r in and r out are respectively the inner and the outer radius of the disk, and k and w 



2. The surface density splitting method 




(1) 
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are respectively defined by 

k ~(a + r) [2) 

and 

ZD = 3 

a + r 

Finally, K (respectively E) is the complete elliptic integral of the first (resp. second) 
kind. Note that the integrand of equation (1) diverges as the modulus k — > 1, namely when 
a — > r. The nature of the singularity is twofold: a logarithmic singularity due to the K- 
function, and a hyperbolic one due to the term 1/w. However, this singularity is integrable. 
There are different ways to estimate improper integrals. Here is a simple and efficient recipe. 
Let us write the surface density as 

S(a) = £(r) + <5£(a,r), (4) 



where £(r) is the surface density at the field point and ££(a, r) is the "residual" surface 
density which, by construction, equals zero when a = r. Thus, equation (1) can be written 
as the sum of two contributions, namely: 



(5) 



where 



.homo. , 



r) 



GE(r) 



E(AQ 
w 



-K(k) 



da 



(6) 



is the radial gravitational field due to a radially homogeneous disk and 
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(7) 



corresponds to the residual profile. The point is that the expression for g^ omo ' also exists 
in a closed analytical form (see Appendix A). This quantity is finite everywhere inside the 
disk, even at the disc edges provided that £(r) vanishes continuously there. It is important to 
mention that this restriction is no more necessary if we consider the disk as a tri-dimensional 
system. Besides, g T r cs ' can be easily computed numerically because both <5£(a, r)K(/c) and 
SH(a,r)E(k)/w are fully regular whatever a and r, and especially for a = r (see Appendix 
B for a proof). As a consequence, the accuracy on g£ es - (and subsequently on the total field 
g(r)) depends only on the performance of the quadrature scheme used to performed the 
radial integral in equation (7). 
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In figure 1, we illustrate this simple and efficient technique by considering a disk with 
inner edge r in = 1, outer edge r out = 100 and surface density S(r) oc exp(— r/r out ) as 
commonly used to model spiral galaxies (Freeman, 1970). The graph shows the residual 
profile <5E(a, r), as well as the two regular functions <5£(a, r)K(/c) and 5H(a,r)'E(k)/w. Here, 
the field point where stands the singularity is arbitrarily set to r = 50 (about the middle of 
the disk). 



3. Rotation curve of the disk 



We can now apply this method to compute the rotation law Q(r) of any disk, whatever 
its size and surface density profile. Here, we simply assume that the contribution to the 
rotation is only gravitation (and neglect pressure and viscosity effects 1 ). Thus, Q(r) is 



defined by the relation f2 2 (r) 
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Using the previous expressions for the radial field, these two terms respectively read: 
4GS(r) fr out 
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and 
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(11) 



These two expressions are easy to compute numerically, for reasons mentioned above. 
Let us remind that at the edges, f2 2 es (r) vanishes because E(r) = as assumed (see above). 



4. Concluding remarks 

We have reported a simple and efficient method to compute the radial field and the 
rotation curve of any razor-thin axi-symmetric disk in the mid-plane, using complete elliptic 



1 According to Sofue & Rubin (2001), this is a good approximation because of the small velocity dispersion 
of interstellar gas. 
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Fig. 1. — Residual surface density profile <5£(a,r), and functions <5£(a, r)K(/c) and 
5Y,(a,r)'E(k)/xz> appearing in Eq.(7) (re-scaled for clarity), for a razor-thin disk extend- 
ing from r in = 1 to r out = 100, and with surface density profile E oc exp(— r/r out ). The field 
point is located at r = 50. 
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integrals. As a matter of fact, a similar treatment can be considered for the gravitational 
potential. Even, this splitting technique can be employed to disks with finite (non-zero) thick- 
ness as suggested in Hure (2003), and to the fully tri-dimensional systems as well (Pierens 
& Hure 2004). 
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A. Radial field due to a a razor-thin disk with constant surface density 



Let us begin with the expression for the radial field due to a disk with constant surface 
density, i.e. equation (6). This expression can be written as follows: 
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We can now set u = - in the first integral and v — - in the second one. Changing the 
modulus in the elliptic integrals according to (Gradshteyn & Ryzhik, 1980): 
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with u' 2 = 1 — u 2 and -u < 1, then equation (Al) becomes: 
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Each integral in this expression is known. Actually, we have (Gradshteyn & Ryzhik, 
1980): 

uE S^Ldu = K(u) - E(u) (A5) 
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Inserting these relations into equation (A4) and re-arranging terms yields the formula 
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B. Singularity removal 



In order to compute numerically equations (7) and (11), we have to check that SH(a, r)K(k) 
and SH(a,r)E(k)/w remain finite when a — > r, which is not obvious at first glance. Since 
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SH(r,r) = by construction, a Taylor expansion of 5T, around a = r yields: 

5Z(a,r) = (a-r)(^\ + O ((a - rf) 

\ / a=r 

= (« - r) (j?) + O ((o - rf) (Bl) 



(a + r)w 



where (^) a _ r is in general finite everywhere. Since E(/c) is bounded, this expansion 
shows that 5Tl(a, r)E(k)/w is finite at the field point (see also figure 1), namely 

ss{r , r) m = {a+r) ^^ ( B2) 

Besides, since K(/c) — > |ln(16/ct7 2 ) as — > f (Abramowitz & Stegun, 1970), we have 
5E(r, r)K(l) = for a = r (see figure 1 again). 



